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SECTION  I 


INTRODUCTION 

During  the  last  three  decades  i  the  research  and  development  in 
finite  element  method  have  grown  from  infancy  to  maturity  and  this 
method  has  revolutionalized  the  methods  of  structural  analysis  and 
design.  Such  development  has  been  well  documented  in  the  texts  by, 
for  example,  Martin  [1],  Zienkiewicz  [2],  and  Gallagher  [3]. 

During  the  last  decade,  the  development  of  calculators  and  desktop 
microcomputers  has  grown  at  a  very  rapid  rate  and  it  becomes  an  obvious 
trend  that  more  simple  daily  structural  designs  will  be  performed  using 
microcomputers.  To  cope  with  such  trend,  works  have  currently  been 
underway  to  convert  some  general  purpose  finite  element  programs 
suitable  for  microcomputers. 

During  the  last  two  decades,  the  research  and  development  of 
laminated  composite  structures  has  grown  at  an  extremely  rapid  pace  and 
it  becomes  an  obvious  trend  that  more  and  more  composite  material  will 
be  used  in  the  design  of  structures  when  the  weight  and  strength  are  of 
primary  consideration.  The  fundamental  development  in  the  mechanics 
of  composite  materials  has  been  documented  by,  for  example,  Tsai  [4] 
and  Jones  [5].  The  basic  theory  of  the  mechanics  of  composite  materials 
particularly  for  laminated  plates,  has  been  widely  used  in  finite 
element  formulations.  Thus,  it  is  common  that  existing  isotropic  and 
homogeneous  finite  elements  also  have  the  capability  of  treating 
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laminate  composite  materials.  A  survey  of  recent  research  in  the 
analysis  of  composite  plates  including  finite  element  methods  can  be 
found,  for  example,  in  [6]  by  Reddy.  An  evaluation  of  finite  element 
software  for  stress  analysis  of  laminated  composites  was  given  in 
Reference  [7]. 

In  view  of  the  trend  of  growing  use  of  microcomputers  in  the 
common  engineering  office,  it  appears  that  there  is  a  definite  need 
for  research  and  development  work  to  tailor  and  simplify  the  basic 
formulations  for  laminated  composite  finite  elements  into  a  form 
suitable  for  programming  using  microcomputers.  As  a  first  step  to 
respond  to  this  need,  in  Chapter  2,  a  simplified  symmetrically 
laminated  beam-type  finite  element  is  developed  and  programmed  for  a 
microcomputer.  The  element  is  assumed  to  have  six  degrees  of  freedom 
at  each  of  the  two  ends:  transverse  deflection  and  slope  due  to  bending 
and  shear,  respectively;  and  a  twisting  angle  and  its  derivative  with 
respect  to  beam  axis. 

This  program  implemented  in  a  microcomputer  is  capable  of  perfor¬ 
ming  stress  analysis  of  symmetrically  laminated  beam  structures  with  a 
single  or  combined  effect  of  bending  moment,  twisting  moment,  and  shear 
deformation,  and  with  arbitrary  loading  and  boundary  conditions.  It  is 
also  capable  of  performing  free  vibration  analysis  without  shear 
deformation.  For  static  analysis,  the  program  has  the  capability  of 
providing  both  numerical  data  and  graphical  plots  of  the  distributions 
of  displacements,  bending  and  twisting  moments,  ply  stresses,  and  the 
portions  contributed  by  shear  deformation.  For  free  vibration  analysis 
the  program  gives  the  natural  frequencies  and  mode  shapes. 


The  simple  homogeneous  and  isotropic  beam  finite  element  formula¬ 
tion  is  essentially  the  same  as  that  traditionally  formulated  by  the 
slope-deflection  equations  for  beams.  Such  an  element  was  extended  by, 
among  others,  McCalley  [8],  Archer  [9],  and  Kapur  [10],  to  include  the 
effects  of  shear  deformation.  Archer's  element  has  two  degrees  of 
freedom  at  each  of  the  two  nodes:  total  transverse  deflection  due  to 
combined  effect  of  bending  and  shearing  deformations  and  its  derivative 
with  respect  to  beam  axis.  On  the  other  hand,  Kapur  developed  a  beam 
element  where  the  bending  and  shearing  deformations  were  considered 
separately  with  separate  associated  degrees  of  freedom.  Thus  the 
element  has  four  degrees  of  freedom  at  each  of  the  two  nodes:  transverse 
deflection  and  its  derivative  with  respect  to  the  beam  axis  due  to 
bending  and  shearing  deformations,  respectively.  For  the  present 
laminated  composite  beam  finite  element,  Kapur's  type  of  approach  was 
used  to  account  for  the  effect  of  shear  deformation.  The  homogeneous 
anisotropic  beam  theory  which  considered  the  coupling  between  bending 
and  torsion  was  given  by  Lekhnitskii  [11].  Such  coupling  was  incor¬ 
porated  in  the  present  formulation.  A  beam  finite  element  formulation 
for  laminated  composite  material  with  a  single  fiber  orientation  and 
shear  deformation  was  given  by  Teh  and  Huang  [12]  for  free  vibration 
analysis.  In  that  element,  six  deyrees  of  freedom  were  assumed  at  each 
node:  total  deflection,  total  slope,  twist  derivative  of  bending  slope 
with  respect  to  beam  axis,  second  derivative  of  bending  slope  with 
respect  to  beam  axis,  angular  displacment  and  angle  of  twist.  This 
paper  differs  from  Reference  [12]  in  that  different  types  of  degrees  of 
freedom  are  assumed  aiming  at  performing  static  and  free  vibration 
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analyses  using  a  microcomputer  in  a  simpler,  more  efficient  and  general 
fashion. 

In  Chapter  3,  a  simple  yet  efficient  formulation  for  a  laminated 
composite  plate  finite  element  was  developed  and  pertinent  efficient 
numerical  algorithms  were  adopted  so  that  the  development  is  suitable 
for  use  in  stand  alone  desktop  microcomputers  for  structural  analysis 
and  design.  Furthermore,  a  previously  formulated  8  d.o.f.  beam  finite 
element  [13]  (see  Figures  1.1  and  1.2)  were  further  evaluated  and 
compared  with  the  more  sophisticated  plate  finite  element. 

In  the  finite  element  formulations,  an  18  d.o.f.  triangular  plate 
element  in  bending  with  anisotropic  symmetrically  laminated  composite 
materials  was  formulated.  The  triangular  plate  finite  element 
described  in  Figure  1.3  has  6  d.o.f. 's  at  each  of  the  3  nodes: 
transverse  deflection,  its  first  derivatives  with  respect  to  K  and  n, 
respectively,  its  second  derivates  with  respect  to  C  and  n,  respective¬ 
ly,  and  its  twist  derivative  with  respect  to  £  and  n,  where  £  and  n 
are  the  local  coordinates.  This  element  was  developed  by  Cowper  et  al. 
[14]  for  isotropic  materials.  It  has  now  been  extended  to  include  the 
effect  of  laminated  composites  for  use  on  the  microcomputer.  This 
element  is  highly  efficient  as  it  meets  the  convergence  criteria  [15] 
and  the  normal  slope  varies  cubically  along  the  interelement  boundaries 
with  full  compatibility.  Furthermore,  the  integrations  of  the  stiff¬ 
ness  coefficients  can  be  expressed  in  a  closed  form.  Both  features 
made  this  element  well-suited  for  microcomputers. 

In  the  development  of  numerical  procedures,  emphasis  has  been 
placed  upon  the  minimization  and  condensation  of  memory  storage  and 


efficiency  of  computation.  The  minimization  of  memory  storage  is 
achieved  by  the  use  of  banded  and  profile  equation  solvers  for  the 
static  case.  For  the  eigenvalue  problem,  a  symmetrical  matrix  eigen¬ 
value  solver  is  used  to  minimize  memory  storage.  For  the  special  case 
of  lumped  mass  matrices,  a  condensation  technique  eliminating  the  zero 
terms  along  the  diagonal  of  the  mass  matrix  is  used  to  reduce  the  total 
number  of  d.o.f.'s  or  equations.  For  efficiency  of  computation,  only 
a  few  of  the  element  stiffness  matrices  need  be  calculated  provided  the 
structure  has  repeating  elements.  Certain  20x20  transformation 
matrices  can  be  inverted  a  priori  without  the  use  of  numerical  routines 
provided  the  elements  are  right  triangles. 

To  demonstrate  and  evaluate  the  present  formulations  and  numerical 
techniques,  a  series  of  static,  free  vibration,  and  buckling  analyses 
of  anisotropic  symmetrically  laminated  plate  problems  have  been 
performed  using  a  microcomputer.  To  further  test  the  previously  formu¬ 
lated  beam  finite  element  [13],  results  of  the  beam  element  has  been 
obtained  and  compared  to  those  of  the  plate  element.  As  will  be 
demonstrated  in  Chapter  3,  the  present  development  is  simple,  efficient, 
and  practical . 

It  has  been  know  that  the  effect  of  shear  deformation  must  be 
included  in  the  finite  element  formulation.  There  are  a  number  of  ways 
to  include  the  effects  of  shear  deformation.  A  popular  approach  is  to 
use  revised  plate  theory  by  Reissner  [16].  However,  in  some  formula¬ 
tions  the  problem  known  as  "shear  locking"  is  encountered  in  the 
inclusion  of  the  effect  of  transverse  shear  deformation.  Some  techni¬ 
ques  such  as  reduced  integration,  penalty  methods,  energy  balancing. 


mixed  formulations  and  others  have  been  proposed  to  alleviate  the  shear 
locking. 

In  Chapter  4,  an  approach  was  presented  in  the  formulation  of  a 
36  d.o.f.  symmetrically  laminated  composite  triangular  plate  element 
which  does  not  present  the  shear  locking  problem.  This  is  an  extended 
version  of  the  18  d.o.f.  isotropic  triangular  thin  plate  element  by 
Cowper,  Kosko,  Lindberg,  and  Olson  [14].  This  formulation  is  also  a 
generalization  of  the  approach  by  Kapur  [10]  for  an  isotropic  beam 
element  with  the  effect  of  shear  deformation.  In  this  approach,  the 
effect  of  shear  deformation  was  incorporated  into  the  formulations  by 
expressing  the  strain  energy  of  deformation  as  the  sum  of  the  thin 
plate  flexural  energy  and  the  energy  of  transverse  shear  deformation, 
and  also  by  expression  the  total  transverse  displacement  as  the  sum  of 
the  displacements  due  to  bending  alone  and  that  due  to  shear  deforma¬ 
tion  alone.  Furthermore,  only  the  displacement  due  to  bending  alone 
appears  in  the  flexural  strain  energy  and  the  displacement  due  to  shear 
deformation  appears  in  the  shear  strain  energy.  A  similar  approach 
for  the  inclusion  of  transverse  shear  deformation  has  been  presented 
by  Bhashyam  and  Gallagher  [17]  for  the  static  analysis  of  isotropic 
plates  using  an  18  d.o.f.  triangular  finite  element.  In  their  approach 
the  strain  energy  expression  is  the  sum  of  the  bending  and  shear  defor¬ 
mation  energies.  The  strain  energy  expression  is  in  terms  of  the  total 
displacement  and  that  due  to  bending  alone.  Only  the  displacement  due 
to  bending  alone  appear  in  the  flexural  strain  energy.  However,  the 
displacement  due  to  shear  deformation  in  the  shear  strain  energy  is 
expressed  as  the  difference  of  the  total  displacement  and  the 
displacement  due  to  bending. 


To  evaluate  the  range  of  applicability  and  efficiency  of  the 
present  formulation,  a  variety  of  isotropic,  orthotropic  and  aniso¬ 
tropic  plates  with  various  thickness-to-length  ratios  were  obtained 
and  compared  with  various  alternative  solutions  available  in  the 
literature  to  demonstrate  this  approach  for  the  free  vibration  of 
plates  with  the  effect  of  shear  deformation  and  rotatory  inertia. 


SECTION  II 


STATIC  AND  DYNAMIC  FORMULATION  OF  A  SYMMETRICALLY  LAMINATED 
BEAM  FINITE  ELEMENT  FOR  A  MICROCOMPUTER 

1.  Symmetrically  Laminated  Finite  Element  Formulation 

In  this  formulation,  a  symmetrically  laminated  composite  beam  was 
considered.  The  beam  is  made  of  layers  of  orthotropic  material  in 
which  the  orthotropic  axes  of  layer  may  be  oriented  at  an  arbitrary 
angle  with  respect  to  the  beam  axis  and  in  which  the  layers  are  stacked 
symmetrically  with  respect  to  the  midsurface  of  the  beam  (see  Figure 
1.1).  For  the  case  of  a  symmetrically  laminated  composite  beam,  there 
exists  coupling  between  the  bending  and  twisting  curvatures. 

1.1.  Formulation 

The  stress-strain  relation  due  to  bending  deformation  only  can  be 
written  as 


where  the  o's  are  the  stress  components,  e's  are  the  strain  components, 
and  the  Q.j's  are  the  inplane  ply  stiffnesses. 
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The  stress-strain  relation  due  to  transverse  shear  deformation  is 


I 


given  by 


°xz  ~  ^44  Exz 


(2.2) 


For  the  general  laminate,  the  linear  strain  distribution  across  the 
thickness  of  the  laminate  is  given  by 
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0  , 

e  +  e 

y  y 
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xy  xy 


(2.3) 


where  e  's  are  the  strains  due  to  inplane  forces,  and  e's  the  strains 
due  to  bending  deformation.  The  moment-curvature  and  inplane  stress- 
strain  relations  for  the  general  laminate  is  given  by 
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where  N's  are  the  inplane  forces,  M’s  the  bending  and  twisting  moments, 


m 


k's  the  bending  and  twisting  curvatures,  A-j's  the  inplane  moduli, 

B. .'s  the  coupling  moduli,  and  D-.'s  the  flexural  moduli. 

'  J  * 

For  the  case  of  a  symmetrical ly  laminated  plate,  the  coupling 
moduli  B,  •  vanish.  Thus  the  flexural  and  membrane  behaviors  become 
uncoupled. 

For  the  case  of  the  beam,  no  bending  moment  exists.  It  is 
noted,  however,  that  the  curvature  k^  is  assumed  to  be  nonzero.  Thus, 
the  moment  curvature  relation  may  be  rewritten  as 

(2.5) 


The  transverse  shear  force-strain  relation  is  given  by 

Qx  =  «D44yx  (2.6) 

where  Qx  is  the  transverse  shear  force,  a  the  shear  coefficient  *  °44 
the  transverse  shear  stiffness  and  yx  the  transverse  shear  strain. 

The  strain  energy  expression  is  given  as 

U  =  I;  f { M  k  +  M  k  }  dA  +  5  [{Q  y  }dA  (2.7) 

2  J  x  x  xy  xy  ?  J  Xx'x  ' 

2 

fi  Wh 

where  k  =  — ,  k  =  2  d$/dx, 

*  dx^  *y 

y  =  dW  /dx,  b  -  width  of  beam,  A  =  area 

A  J 
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=  displacement  due  to  bending  deformation. 


W$  =  displacement  due  to  shear  deformation, 

4  =  twist. 

The  kinetic  energy  expression  for  a  plate  can  be  written  in  the 
general  terms  as 


T  =  ^ 
1  2 


U 

4  4 


W  (x,y,t)dxdydz 


(2.8) 


where  p  is  the  mass  density  per  unit  volume.  The  dot  represents 
derivative  with  respect  to  time. 

For  the  case  of  the  beam  in  the  absence  of  shear  deformation,  the 
deflection  function  may  be  written  as 


(2.9) 

! 

I 


W(x,y)  =  Wb(x)  +  y  <j>(x) 


Substituting  Equation  (2.9)  into  Equation  (2.8)  and  performing  integra¬ 
tion  result  in. 


I 


£ 

f  W2  (x)  dx  +  ^  ;2(x)  dx 

'  n 


(2.10) 


where  A  =  cross-sectional  area  of  the  beam,  and  J  =  polar  mass  moment 
of  inertia  about  the  centroidal  axis. 

Substituting  the  kinetic  energy  expression  Equation  (2.10)  and  the 
strain  energy  expression  Equation  (2.7)  into  the  Lagrange's  equation 
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j-W/.wi  ...i-i 


{F}  =  [k]{q}  +  [m] { q } 


(2.12) 


where  [k]  is  the  stiffness  matrix;  [m]  is  the  consistent  mass  matrix, 
(which  can  also  be  derived  in  other  forms  such  as  "lumped");  and  {q} 
is  the  vector  for  the  six  nodal  degrees  of  freedom. 

Assuming  free  vibration  with  simple  harmonic  frequency  u,  Equation 
(2.12)  becomes 

{[k]  -  iu2[m]}  {q}  =  0  (2.13) 


The  consistent  mass  matrix  and  the  stiffness  matrix  including  shear 
deformation  are  given  by 
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where 


'i'l 

♦M 


111  =  156  F;  M12  =  -22L  F;  M15  =  54  F;  M16  =  13L  F;  H22  =  4L"  F; 

125  -  -13L  F;  M26  =  -3L3  F;  H33  =  156  G;  H34  =  -22L  G;  M37  =  54  G; 

<38  =  13L  G;  M55  =  156  F;  M56  -  22L  D;  M66  =  4L2  F;  M77  =  156  G; 

178  =  22L  G;  M88  =  4L2  G;  G  =  JL/420;  F  =  mL/420; 

)  *  mass  polar  moment  of  inertia  per  unit  length  of  the  element; 
n  =  mass  per  unit  length  of  the  element,  and 
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where 

Kll  =  - K 1 7  =  K77  =  12  D11*L3;  K12  =  K18  =  -K27  =  -K78  =  -6  D1 1*/L2 
K22  =  K88  =  4  Dll*/L;  K28  =  2D11*/L;  K16  =  -K1C  =  -K25  -  K2B  =  K56 
-K67  =  K7C  =  -K8B  =  6  066V5L;  K56  =  K5C  =  -KBC  =  -D66*/10;  K66  = 

2  066*  1/15;  K6C  -  066*  L/30;  K33  *  -K39  =  K99  =  6  S/5L;  K34  =  K37 
-K9A  =  -K49  =  -A/10;  K44  =  KAA  =  2L  S/15;  K4A  =  -L  S/20;  S  =  aD44; 


Dll*  =  b  (D11-D12VD22);  D16*  =  b(D16-D12  D26/D22);  D66*  =  b(D66-D26V 
D22). 


1.2  Description  of  the  Element 

The  present  symmetrically  laminated  beam  element  is  described  in 
Figure  2.1.  The  element  possesses  6  degrees  of  freedom  at  each  of  the 
two  nodes:  the  deflection  due  to  bending  Wb>  the  deflection  due  to 
shear  deformation  W$,  and  their  respective  derivatives  with  respect  to 
the  x-axis  eb  (=  -dW^/dx)  and  9S  (=  -dWs/dx),  and  the  twisting  angle 
4>  (=t)  and  its  derivative  t ‘ ( =d4>/dx ) .  The  displacement  functions  for 
W^,  and  4>  are  assumed  as, 

wb(x)  =  f1(x)wbi  +  f2(x)ebi  +  f3(x)wb2  +  f4(x)ob2 

W$(x)  =  FjMW  +  f2(x)e  +  f3^x^Ws2  +  f4^x^6s2  (2.16) 

♦(x)  =  f1(x>T1  +  f2(x)Ti  +  f3(x)T2  +  f4(x)x2 

where  the  shape  functions  are  in  terms  of 

fj(x)  =  l+2(x/£)3=3(x/£)2, 
fp(x)  =  -x+2x2/£-x3/£2 , 

(2.17) 

f 3 ( x )  =  3(x/£)2-2(x/£)3, 
f4(x)  -  -x3/£2+x2/£. 
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2.  Microcomputer  Program 

The  formulation  of  the  present  12  d.o.f.  symmetrically  laminated 
composite  beam  finite  element  has  been  coded  into  a  microcomputer  in 
Basic  language.  For  static  analysis,  the  program  computes  the 
transverse  deflections  and  its  slope  distribution  due  to  bending  and 
shearing  deformation  along  the  beam  axis.  It  then  computes  the  moment 
shear  force,  and  torque  distributions  along  the  beam  axis.  It  finally 
computes  the  inplane  normal  and  shearing  stresses,  and  transverse 
shearing  stress  distributions  through  the  laminate  thickness.  For 
free  vibration  analysis,  the  program  computes  the  natural  frequencies 
and  the  corresponding  mode  shapes  along  the  beam  axis.  A  flowchart  of 
the  program  is  given  in  Figure  2.2. 

For  the  static  case,  the  program  uses  a  symmetrically  banded 
matrix  solver  which  reduces  both  computing  time  and  memory  storage. 

The  program  also  has  a  graphics  routine  which  plots  the  distributions 
of  the  various  aforementioned  quantities  on  the  microcomputer  monitor 
screen.  There  are  several  devices  which  can  be  used  to  accelerate  the 
computational  speed.  The  first  method  is  to  compile  the  program  using 
the  available  compilers  on  the  market.  The  second  method  is  to  use  a 
machine-dependent  arithmetic  chip  which  accelerates  to  store  the 
results  in  a  file,  which  can  be  printed  later.  Similarly,  the  plots 
may  be  stored  in  a  file  and  recalled,  and  a  hard  copy  of  the  plot  can 
be  obtained  using  available  "screen  dumps." 

For  the  free  vibration,  the  program  uses  a  Jacobi  eigenproblem 
solver.  It  computes  the  natural  frequencies  and  plots  the  mode  shapes 


USER  INPUT: 

GEOMETRY,  MATERIAL  PROPERTIES,  AND  LAMINATE  LAYUP 


PRINTED  OR  GRAPHIC  OUTPUT: 

STATIC  CASE:  DEFLECTIONS,  TWISTING  ANGLES,  SHEAR  FORCE, 
BENDING  MOMENT,  TWISTING  MOMENT  ALONG  THE 
BEAM,  PLY  STRESSES  THROUGH  THE  THICKNESS, 
AND  THE  PORTIONS  CONTRIBUTED  BY  SHEAR 
FREE  VIBRATION  CASE:  NATURAL  FREQUENCIES  AND  ITS 

CORRESPONDING  MODESHAPES 


Figure  2.2  Flow  chart  for  the  microcomputer  program. 


This  program  has  been  written  for  the  IBM  PC.  The  minimum  config¬ 


uration  needed  to  run  the  program  requires  a  128K  memory  storage,  a 
single-sided  disk  drive,  and  a  monitor.  Additional  hardware  and 
software  can  enhance  the  performance.  The  compiler  provided  by  IBM  is 
the  IBM  Basic  Compiler  Version  1.00.  The  additional  hardware  is  the 
8087  arithmetic  chip  which  is  used  in  conjunction  with  the  compiler. 

The  "screen  dump"  software  is  provided  by  the  DOS  System  Diskette. 

3.  Evaluative  Analysis 

The  formulation,  solution  procedure,  and  program  have  been 
evaluated  by  performing  the  following  examples  with  the  alternative 
solutions  for  comparison. 

3.1.  Static  Analysis 

3.1.1.  Homogeneous,  Isotropic  Cantilever  Beam  Linder  End  Load  P 
To  test  the  portion  of  the  formulation  which  accounts  for  the 

effect  of  shear  deformation,  an  example  of  a  homogeneous,  isotropic 
cantilever  beam  under  end  load  P  was  first  analyzed.  The  results 
obtained  using  one  12  d.o.f.  element  for  the  end  deflection  due  to  the 
effect  of  shear  deformation  only  are  given  in  Table  3.1  for  two 
different  values  of  shear  coefficient  a,  as  defined  in  Reference  [18]. 
For  a  =  0.667,  an  alternative  exact  solution  available  in  the  text  by 
Timoshenko  [19]  and  a  one-element  solution  by  Archer  [9]  are  shown  in 
the  table.  The  four  values  for  the  end  deflection  by  the  four  different 
methods  are  seen  to  be  in  excellent  agreement.  For  a  =  0.867,  a 
one-element  by  Archer  [13]  and  a  solution  based  on  the  energy  method  by 
Popov  [20]  are  also  shown  in  the  table;  the  present  solution  is  4 % 


Table  2.1  Maximum  deflection  W$  in  Pch  /4EI  due  to  shear 

deformation  only  for  a  homogeneous  isotropic 
cantilever  beam  under  end  load. 


lower  than  that  by  Popov  and  17%  lower  than  that  by  Archer.  It  is 
noted  that  the  difference  between  the  present  solution  and  that  by 
Archer  may  be  due  to  the  difference  in  the  boundary  conditions  at  the 
fixed  end.  In  the  former,  it  is  assumed  that  dW$/dx  /  0,  and  in  the 
latter,  it  is  assumed  that  dW/dx  +  <i>  =  0,  where  W$,  is  the  deflection 
due  to  shear  deformation;  W  is  the  total  deflection;  and  ^  is  the  shear 
angle. 

3.1.2.  Simply-Supported  Rectangular  3  Layer  ( 0/90/0 )  Cross-Ply 
Laminated  Plate  With  Infinite  Aspect  Ratio  (Wide  Beam)  With  Shear 
Deformation 

The  second  example  chosen  was  a  simply  supported  rectangular  3 
layer  (0/90/0)  cross  ply  laminated  plate  with  infinite  aspect  ratio 
(wide  beam)  subjected  to  a  sinusoidally  distributed  load, 

P  =  PQ  sin  (2.18) 

where  t  is  the  length  and  x  is  the  coordinate  along  the  span.  Due  to 
symmetry,  only  half  of  the  beam  need  be  modeled;  and  1,2, 3, 6,  and  8 
elements  were  used,  respectively.  The  work  equivalent  loads  based  on 
the  integration  of  the  products  of  the  shape  functions  and  the  distri¬ 
buted  load  were  used.  The  central  deflection  due  to  the  effect  of  shear 
deformation  only  was  obtained  and  given  in  Table  2.2 

An  alternative  analytical  solution  provided  by  Pagano  and  Whitney 
[21]  is  also  shown  in  Table  2.2  It  is  seen  that  the  present  solution 
practically  converged  to  that  given  in  Reference  [2l]  at  the  6  element 
level . 
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Table  2.2  Non-dimensional ized  maximum  deflection  due  to  shear 
deformation  only  of  a  simply-supported  rectangular 
3  layer  (0/90/0)  cross-ply  laminated  plate 
with  infinite  aspect  ratio  (wide  beam)  under  a 


sinusoidal  load  P 


sin 
o 


^ ) 
.«Ki 


Number  of  Elements 
for  half  of  beam 


Non-dimensional ized  max¬ 
imum  deflection 

Whitney  (21) 

V  D44 

\  °44 

v2 

v2 

0.3599 

0.4463 

0.4014 

0.4053 

0.4057 

i 

0.4056 

R 


3.1.3.  Anisotropic  16  Layer  (45,,/-45,,  )s  Laminated  Cantilever  Beam 


Under  End  Load  P 

The  third  example  chosen  was  an  anisotropic  16  layer  laminated 
cantilever  beam  under  an  end  load  P  without  the  effect  of  shear 
deformation.  The  results  for  the  distributions  of  twisting  angle  and 
deflection  due  to  bending,  respectively,  along  the  beam  length  were 
obtained  used  one  element  only  and  were  plotted  in  Figure  2.3.  With 
the  use  of  the  displacement  functions  in  Equation  (1),  the  displace¬ 
ment  can  be  interpolated  anywhere  within  the  element.  The  twisting 

o 

angle  and  the  deflection  were  nondimensional ized  as  (4t/P i  djg)  and 

(3W/P£  d^),  respectively,  where  t  =  twisting  angle;  l  =  length  of  the 

beam,  d^,  d^  are  flexural  compliances;  and  W  =  deflection. 

An  alternative  analytical  solution  for  this  problem  treating  the 

beam  as  homogeneous  and  anisotropic  by  Lekhnitskii  [11]  is  also 

plotted  in  Figure  2.4  for  comparison.  The  agreement  is  good. 

The  normalized  inplane  ply  stress  (o  ,  o  ,  and  o  )  distribu- 

xx  xy  yy 

tions  through  the  laminate  thickness  were  plotted  in  Figure  2.4  It  is 

noted  that  force  equilibrium  is  satisfied,  and  the  bending  moment  due 

to  end  load  P  is  recovered  from  summing  the  moments  due  to  o  .  It  is 

also  noted  that  the  summation  of  the  moments  due  to  o  and  o  , 

xy  yy’ 

respectively,  give  zero  resultant  moments. 


3.1.4.  Quasi-Isotropic  (0/90/‘45)s  Laminated  Beam  Under  Four  Point 


The  fourth  example  chosen  was  a  quasi-isotropic  ( 0/ 90/ 4 4 5 ) s 
laminated  beam  under  four  point  bending.  The  results  for  the  distribu¬ 
tion  of  the  transverse  shear  stress  through  the  laminate  thickness  were 


NON-DIMENSIONAL  BEAM  DISPLACEMENT 


BEAM  LENGTH  j 


Figure  2.3  Distribution  of  deflection  and  twisting  angle 
due  to  bending  for  a  16  ply 
anisotropic  laminated  cantilever  beam  under 
end  load  P. 
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end  load 


obtained  using  two  elements,  and  were  plotted  in  Figure  2.5  The  shear 
stress  has  been  non-dimensional ized  as  (°x2l/Qxh  )*  The  resultant 
shear  force  calculated  from  under  the  shear  stres'-  curve  is  1.03,  3% 
higher  than  unity.  The  result  given  by  Whitney  [22]  is  also  shown  in 
the  figure. 

3.1.5.  Anisotropic  16  Layer  ( 45^/-454 ) s  Laminated  Cantilever  Beam 
Under  End  Load  P  With  Shear  Deformation 

The  final  example  chosen  was  the  same  as  that  given  in  (2.3. 1.3) 
but  with  the  effect  of  shear  deformation.  The  results  for  the  distri¬ 
butions  of  twisting  angle  and  deflection  due  to  bending  as  well  as 
shear  deformation  along  the  length  of  the  beam  were  obtained  using  one 
element  and  were  plotted  in  Figure  2.6.  Again,  the  twisting  angle  and 
the  deflections  were  nondimensinal ized  the  same  way  as  were  those  done 
in  Figure  2.3.  The  ply  stresses  are  the  same  as  those  in  Figure  2.4 

3.2.  Free  Vibration  Analysis 

3.2.1.  Homogeneous,  Isotropic,  Thin  Cantilever  Plate  Without  Shear 
Deformation 

To  test  the  portion  of  the  formulation  which  accounts  for  the 
special  case  of  the  homogeneous,  isotropic  materials,  an  example  of 
an  aluminum  cantilever  plate  without  shear  deformation  was  first 
analyzed.  The  results  obtained  using  four  present  elements  are  given 
in  Table  2.3.  The  experimental  results  used  for  comparison  are  given 
by  Crawley  [23].  The  present  shows  good  agreement  for  the  bending 
frequencies  is  due  to  the  present  modeling  of  the  two-dimensional  plate 


as  a  one-dimensional  beam. 
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Figure  2.6  Distribution  of  deflections  due  to  bending  and 
shear  deformation  and  twisting  angle  for  a  16 
ply  (45^/-45^)s  anisotropic  laminated  cantilever 

beam  under  end  load  P. 
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Table  2.3  Natural  frequencies  of  an  isotropic,  homogeneous 
aluminum  cantilever  plate  (wide  beam)  without 
shear  deformation. 


3.2.2.  Thin  Anisotropic  6  Layer  [(^/Ojs  Laminated  Cantilever  Plates 
Without  Shear  Deformation 

An  anisotropic  laminated  plate  with  [(^/Ojs  stacking  sequences  was 
analyzed.  The  results  obtained  using  a  variety  of  numbers  of  the 
present  element  are  given  in  Table  2.4,  for  both  consistent  and  lumped 
mass  formulations.  The  experimental  results  used  for  comparison  are 
given  by  Jensen  [24].  The  present  solution  is  seen  to  converge  at  the 
four  element  level  and  to  wel 1 -converge  at  the  eight  element  level. 

The  present  solution  is  in  good  agreement  with  the  experimental 
results  for  the  bending  frequencies.  The  percentage  of  discrepancies 
between  the  presently  obtained  torsional  frequencies  and  the  experimen¬ 
tal  values  are  within  the  approximate  range  of  10  -  20%.  Again,  this 
discrepancy  may  be  due  to  the  present  modeling  of  the  two-dimensional 
plate  as  a  one-dimensional  beam. 
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ng,  second  bending,  and  first  torsional  modes. 


SECTION  III 


r 

t 

’« 

4 


» 

t] 

I 

4 


STATIC,  DYNAMIC,  AND  BUCKLING  FORMULATION  OF  A  SYMMETRICALLY 
LAMINATED  PLATE  FINITE  ELEMENT  FOR  A  MICROCOMPUTER 

1.  Formulations 

The  beam  conventions  are  given  in  Figure  1.1.  The  8  d.o.f.  beam 
finite  element  formulated  previously  is  given  in  Figure  1.2. 

The  plate  element  described  in  Figure  1.3  was  developed  for  the 
isotropic  case  by  Cowper  et  al.  The  details  of  the  derivations  can  be 
found  in  [15,  25]  so  they  will  not  be  presented  here.  Only  those 
pertaining  to  the  extension  of  the  element  to  include  the  properties 
of  composite  materials  are  given  here. 

The  strain  energy  expression  can  be  written  as  follows  for  a 
symmetrical ly  laminated  composite  plate 


where  Wrf,  W  ,  2Wr  are  the  bending  and  twisting  curvatures  in  local 
coordinates  £  and  n  of  the  elements,  and  D.j  are  the  flexural  moduli. 
The  integration  is  over  the  area  of  the  triangular  element. 
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2.  Microcomputer  Program 

The  formulations  of  the  present  18  d.o.f.  symmetrically  laminated 
composite  plate  finite  element  has  been  coded  into  a  microcomputer  in 
FORTRAN  language.  The  advantage  of  using  FORTRAN  language  over  BASIC 
are  two-fold.  The  first  is  that  the  program  can  be  developed  more 
rapidly  on  a  mainframe  computer  and  then  downloaded  to  the  micro¬ 
computer.  The  second  is  that  many  algorithms  already  developed  are  in 
FORTRAN  and  these  can  be  readily  adapted  and  modified  for  use  on  the 
microcomputer.  There  now  exist  many  commercial  FORTRAN  compilers  for 
microcomputers. 

Many  algorithms  already  available  have  been  adopted  and  modified 
for  use  on  the  microcomputer.  Those  algorithms  that  have  been  used  in 
the  plate  program  are  presented  here. 

For  static  analysis,  a  banded  [26]  and  a  skyline  [27]  equation 
solver  is  used  to  minimize  memory  storage.  The  boundaries  of  the 
banded  and  skyline  storage  are  illustrated  in  Figure  3.1.  As  an 
example,  a  comparison  of  various  storage  scheme  for  the  stiffness 
matrix  for  a  simply  supported  square  isotropic  plate  under  uniform  load 
is  given  for  3  different  meshes  in  Table  3.1.  For  eigenvalue  analysis, 
a  symmetric  matrix  eigenvalue  solver  is  used  to  minimize  memory 
storage.  The  Jacobi  eigenvalue  solver  [28]  has  been  modified  to  use  a 
symmetric  matrix  storage  scheme. 

There  are  two  methods  of  lumping  masses.  The  first  method,  let  it 
be  called  lumped  mass  method  1,  can  be  found  in  Reference  [29].  The 
total  mass  of  the  element  is  distributed  in  accordance  with  the  distri¬ 
bution  of  the  mass  along  the  diagonal  of  the  consistent  mass  matrix. 

IS 


Table  3.1  Various  schemes  for  memory  storage  for  [ K] 
for  a  simply  supported  square  isotropic 
plate  under  uniform  load. 


Mesh  for 
a  Quadrant 

Degrees  of 
Freedom 

Storage  Modes 

Full  Matrix  Symmetric  Band  Skyline 

3x3 

4x4 

5x5 

60 

104 

160 

3600  1830  1395  1001 

10816  5460  3569  2737 

25600  12880  6664  5257 

The  second  method  is  based  on  the  most  common  method  of  lumping;  let 
it  be  called  lumped  mass  method  2.  The  total  mass  of  the  element  is 


divided  by  3  and  the  one-third  masses  becomes  the  terms  m^,  m^,  and 
m13-13’  resPect^ve^*  F°r  the  special  case  of  lumped  method  two,  a 
static  condensation  technique  [30]  eliminating  the  zero  terms  along 
the  diagonal  of  the  mass  matrix  is  used  to  reduce  the  total  number  of 
the  equations.  The  eigenvalue  problem  is  given  by 
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(3.2) 


where  ^1  denote  the  degrees  of  freedom  associated  with  mass,  {yQ} 
denote  the  massless  d.o.f.'s,  and  [M]  is  the  diagonal  mass  matrix  with 
non-zero  diagonal  terms. 

The  massless  d.o.f.'s  {yQ}  can  be  expressed  in  terms  of  those 
associated  with  mass  {yp}  as  follows. 

tKp/  {v  = '  CKo°]  <V  <3-3> 

Thus , 


-U2[M](yp>  *  l[Kpp]  -  [Kp0][K00]-'  [Kpo]Tl  lyp)  =  0  (3.4) 


or 


i 

'  38 


(3.5) 


-a»2[M]{yp}  +  [Kn]lyp)  -  0 

where  [K^]  is  a  symmetrical  reduced  matrix.  A  technique  making  use 
of  the  diagonal  matrix  [31]  which  preserves  the  symmetry  of  the 
stiffness  matrix  and  which  allows  use  of  a  symmetric  matrix  eigenvalue 
solver  is  given  below. 

Given 


(3.6) 


and 


[M] 


-  1/2 


l//a 


l//b 


l//c 


l//d 


l//e 


(3.7) 


where  [M]  is  a  diagonal  matrix  with  nonzero  diagonal  terms,  and  [M]~^2 
is  a  matrix  with  inverse  square  root  diagonal  terms. 


then 


[M]’1/2  [M]  [M]' 1/2  =  [I] 


(3.8) 


and 


V) 


-A[l]ly)  ♦  [MJ 


L KTT ]  [M]  (y)  =  0 


(3.9) 


where  [I]  is  the  identity  matrix.  The  stiffness  matrix  [ ]  premul- 
tiplied  and  postmultipl ied  by  [M]  results  in  a  symmetric  matrix. 

For  static  analysis,  the  program  computes  the  various  displace¬ 
ments  (transverse  deflection,  slopes,  and  bending  and  twisting 
curvatures)  at  the  nodes.  For  eigenvalue  analysis,  the  program 
computes  the  eigenvalues  and  the  corresponding  mode  shapes.  The 
program  gives  a  three-dimensional  view  of  the  plate  static  deflection 
and  plate  vibration  and  buckling  mode  shapes. 


3.  EVALUATIVE  ANALYSIS 

The  formulation,  solution  procedure,  and  program  have  been 
evaluated  by  performing  the  analysis  of  a  series  of  example  problems. 


3.1.  Static  Analysis 
Uniformly  Loaded  Cantilever  Beam 

In  the  first  example,  comparison  is  made  for  the  centerline 
deflections  and  twist  of  two  modelings  for  a  uniformly  loaded  canti¬ 
lever  beam  of  aspect  ratio  of  ten  for  various  fiber  angles.  This 
comparison  is  made  to  further  evaluate  the  previous  beam  element  as  to 
its  accuracy  and  efficiency  in  various  applications.  Beam  elements 
can  be  used  to  save  computational  time  if  its  accuracy  are  close  to 
those  of  the  plate,  but  if  deformation  modes  are  required  then  plate 
elements  have  the  advantage.  The  beam  finite  element  has  4  d.o.f.'s 
at  each  of  its  two  nodes:  and  0^,  the  transverse  deflection  and 

its  slope  due  to  bending,  and  t  and  r‘,  the  twisting  angle  and  the  rate 


*\v; v„w.\y 


of  twisting  angle.  As  shown  in  Table  3.2,  the  plate  model  uses  40 
elements  and  the  beam  model  uses  5  elements.  The  dimensions  of  the 
beam  are  10m  in  length  by  lm  in  width  by  0.0026m  in  thickness.  The 


material  properties  of  the  ply  are  E  =  181  GP  .  E  =  10.3  GP  , 

x  d  y  a 

v  =  0.28,  and  E  =  7.17  GP  .  THe  centerline  deflection  at  the  tip 
xy  b  d 

and  at  60%  of  the  length  from  the  fixed  end  are  given  in  Figure  3.9. 

The  results  show  a  maximum  discrepancy  of  15%  between  the  two  modelings 
at  a  fiber  angle  of  30  degrees.  For  the  centerline  twist,  the  results 
are  given  in  Table  3.3.  These  results  show  a  maximum  discrepancy  of 
15%  at  a  fiber  angle  of  30  degrees.  For  the  fiber  angle  of  30  degrees 
in  this  example,  the  coupling  effect  between  bending  and  twisting  is 
most  pronounced  as  predicted  by  the  shear  coupling.  It  appears  that 
the  beam  element  does  not  predict  the  shear  coupling  well. 

3.2.  Vibration  Analysis 

3.2.1.  Thin  Anisotropic  6  Layer  [e2/0]s  Cantilever  Plate 

The  second  example  gives  the  natural  frequencies  of  a  thin 
anisotropic  cantilever  plate  of  aspect  ratio  of  four.  The  dimensions 
and  the  material  properties  of  the  plate  are  given  in  Reference  [32]. 
The  results  given  in  Table  3.4  using  both  the  lumped  mass  and 
consistent  mass  formulation  of  the  triangular  plate  elements  are 
compared  with  those  results  given  in  Reference  [32].  The  first  two 
columns  are  based  on  lumped  mass  method  1.  The  third  and  fourth 
columns  are  lumped  mass  method  2.  The  results  show  that  the  lumped 
mass  method  1  does  not  seem  to  give  improved  results  over  those  of 
lumped  method  2.  Furthermore,  the  three  cases  in  increasing  computa¬ 
tional  time  and  effort  required  are  lumped  mass  method  2,  lumped  mass 
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method  1,  and  consistent  mass  method.  The  lumped  mass  method  2  uses 
fewer  degrees  of  freedom.  Thus  one  can  take  advantage  of  the  condensa¬ 
tion  of  the  mass  matrices  as  described  in  Equations  3. 6-3. 9.  The 
results  using  consistent  mass  formulation  are  given  in  the  third  and 
fourth  columns.  All  the  results  compare  very  well  with  the  experimen¬ 
tal  and  finite  element  (8  node,  40  d.o.f.,  quadrilateral,  mixed, 
hybrid  element)  results  given  in  Reference  [32]. 

3.2.2.  Anisotropic  6  Layer  [e^/O^s  Cantilever  Beam 

The  third  example  given  in  Table  3.5  is  a  comparison  of  two 
modelings  of  an  anisotropic  cantilever  beam  with  an  aspect  ratio  of 
ten.  The  material  properties  are  given  in  Reference  [32].  The 
dimensions  are  10m  in  length,  lm  in  width,  and  0.804mm  in  thickness. 

The  plate  modeling  uses  both  the  20  (1x10  mesh)  and  40  (2x10  mesh) 
elements,  respectively,  while  the  beam  modeling  uses  4  and  10  elements 
respectively. 

It  is  seen  that  the  three  sets  of  results  obtained  by  using  the 
three  formulations  (1.  Beam  element  with  consistent  mass,  2.  plate 
element  with  lumped  mass,  3.  plate  element  with  consistent  mass),  each 
of  which  using  two  selected  meshes,  are  in  good  agreement.  The 
maximum  discrepancies  occur  in  the  case  with  fiber  angle  of  30  degrees 
where  the  bending  and  twisting  moments  are  more  pronouncedly  coupled. 

In  that  case,  the  discrepancy  in  natural  frequency  of  the  third  mode 
between  the  10  elements  and  40  plate  elements  (both  with  consistent 
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3.3.  Buck! ing  Analysis 

Simply  Supported  Anisotropic  Square  Plate  Under  Compressive  Axial  Load. 

The  dimensions  and  stiffness  properties  of  the  present  plate 
example  are  given  in  Reference  [33].  The  critical  buckling  loads 
obtained  for  various  fiber  orientations  using  32  triangular  elements 
(4x4  rectangular  mesh)  are  shown  in  Table  3.6  The  results  obtained 
by  using  the  Galerkin  method  [33],  and  the  experimental  results  [34], 
and  the  exact  solutions  for  the  orthotropic  case  only  [35]  are  all 
shown  for  comparison.  It  is  seen  that  for  the  cases  with  fiber  angles 
of  0  and  90  degrees,  the  present  results  agree  well  with  the  exact 
solution  [35].  For  the  case  with  fiber  angles  of  30,  45,  and  60 
degrees,  the  present  results  are  in  fair  agreement  with  those  given  by 


Table  3.6  Critical  buckling  loads  for  a  simply  supported 
anisotropic  square  plate  under  uniform 
compressive  axial  load. 


Present 
32  el 
94  dof 


0 

318.86 

285 

271 

318.91 

30 

387.90 

425 

399 

-- 

45 

354.67 

406 

364 

— 

60 

355.40 

381 

433 

— 

90 

203.88 

210 

251 

203.75 

Ref.  [33]  Exp.  [34] 


Exact  [35] 


SECTION  IV 


A  36  DOF  SYMMETRICALLY  LAMINATED  TRIANGULAR  ELEMENT 
WITH  SHEAR  DEFORMATION  AND  ROTATORY  INERTIA 

1.  Formulation 

The  effect  of  transverse  shear  deformation  is  incorporated  into 
the  formulations  by  expressing  the  strain  energy  of  deformation  as  the 
sum  of  the  thin  plate  flexural  energy  and  the  energy  of  transverse 
shear  deformation. 


U  =  U.+  U  (4.1) 

b  s 

The  effects  of  rotatory  inertia  is  included  in  the  formulation  of  the 
kinetic  energy  expression. 

FGr  the  present  plate  finite  element,  two  displacement  functions 
are  used  which  include  the  displacements  due  to  bending  and  that  due 
to  transverse  shear  deformation,  W$.  The  total  displacement  is 

W  =  W.  +  W  (4.2) 

tbs 

Strain  energy  car,  be  expressed  as  follows 


49 


m 


where  the  curvatures  are  given  as 


kx  =  ^V3*2,  ky  =  a2wb/ay2*  and  kxy  =  232Wb/3x3y  (4.4) 


£ 


and  the  shearing  angles  are  given  as 


ux  =  3Ws/Dx  and  uy  =  3W$/3y 


(4.5) 


The  matrices  [D]  and  [A]  contains  the  flexural  and  shearing  moduli, 
respectively. 

It  is  seen  that  the  flexural  energy  is  only  in  terms  of  the 
displacement  due  to  bending  and  that  the  shear  strain  energy  is  only 
in  terms  of  the  displacement  due  to  shear  deformation  W$. 

The  kinetic  energy  including  the  effect  of  rotatory  inertia  is 
given  by 


T  =  p[jj  W2(x,y)dxdy 


+  \  ff  [Wk  +  «k  Jdxdy] 


(4.6) 


vWKv; 


WWi.'W 


where  p  is  the  density  per  unit  volume,  is  the  total  deflection, 

W.  and  W.  are  the  first  derivative  of  the  deflection  due  to  bending 
d,x  n,y 

alone,  and  h,  and  I  are  thickness,  and  moment  of  inertia,  respectively. 


The  displacement  functions  for  the  present  36  d.o.f.  triangular 


element  (see  Figure  4.1)  as  described  in  local  coordinates  are 


20  m.  o,  . 

W.  =  l  a-£  n  m. ,  n.  =  0  for  i  >  20) 
d  i=i  1  11 


(4.7) 


,ri  Si 


W  =  l  a • £  n  r.,  S-  =  0  for  i  <  21 
s  i  =-21  1  1  1 


where  the  displacement  due  to  bending  is  described  in  Reference  [25]. 


The  displacment  due  to  shear  deformation  Ws  is  assumed  to  be  of  the 


same  form  as  that  for  Wb.  Thus  these  displacement  functions  are 


independent  of  each  other. 


The  derivation  of  the  stiffness  matrix  due  to  bending  alone  has 


been  derived  previously  for  the  isotropic  case  by  Cowper  et  al.  [25] 


and  extended  to  the  anisotropic  case  by  the  present  authors  [36]  for 


the  case  of  thin  plate.  All  the  explanations  and  derivations  given  in 


[25]  and  [26]  will  not  be  repeated  here.  Only  the  derivation  of  the 


stiffness  matrix  due  to  shear  deformation  alone  is  given  here. 


The  strain  energy  expression  in  local  coordinates  due  to  shear 


deformation  can  be  expressed  as 


s  "  2 


A44  0 


0  Ar 


dr.dn  (4.8) 


where  W  and  W  are  the  shear  angles  in  local  coordinates  f,  and  n  of 
£  n 

the  elements,  and  A.  are  the  shear  moduli.  Ihe  integration  is  over 

the  area  of  the  triangular  element.  The  expression  in  terms  of 

global  coordinates  follows  the  same  procedure  for  transformation  and 

rotation  of  the  local  coordinates  as  described  in  [2]. 

Shear  correction  factor  and  02  are  incorporated  into  the  shear 

moduli  as  follows  A„.  =  a.  G  h,  and  Acc  =  a0  G  h,  where  G  =  G 

i  yz  dd  c  yz  yz  xz 

is  assumed  for  the  shear  stiffnesses.  They  are  assumed  to  be  5/6 
unless  otherwise  specified.  This  value  was  used  by  Reissner  [16],  and 
Whitney  [37]. 

The  equations  of  motion  for  the  case  of  free  vibration  for  the 
element  can  be  expressed  in  matrix  form  as 


M  M  W. 

s  b 


M  M  W 


0  WL 


0  K_ 


(4.9) 


where 


t«bbi 


consistent  mass  and  matrix  including  the  effect  of 
rotatory  inertia  matrix, 
consistent  mass  matrix, 
bending  stiffness  matrices, 


=  shear  stiffness  matrices. 


degrees  of  freedom  corresponding  to  the  displacement  due 
to  bending,  and 

degrees  of  freedom  corresponding  to  the  displacement  due 
to  shear  deformation ,  and 


oj  =  natural  frequency. 

2.  Evaluative  Analysis 

In  order  to  evaluate  the  performance  of  the  present  finite  element 
formulation,  and  to  find  an  indication  of  the  range  of  applicability,  a 
series  of  symmetrically  chosen  examples  was  analyzed  and  compared  with 
various  alternative  solutions. 

2.1  Isotropic  Simply-Supported  Square  Thick  Plates 

The  first  example  chosen  was  to  find  several  lower  mode  natural 
frequencies  for  an  isotropic  simply  supported  square  plate.  For  this 
example,  an  exact  30  elasticity  solution  [38]  as  well  as  solution  by 
Mindlin  theory  [39]  are  available  for  comparison.  The  results  for  the 
nondimensional  lower  mode  natural  frequencies  for  various  thickness-to 
length  ratio  (h/a)  obtained  using  a  5x5  mesh  of  the  present  36  d.o.f. 
element  are  given  in  Table  4.1.  The  natural  frequencies  were  nondimen 
sionalized  as 


W  =  pi/  3^/ti^D  (4.10) 

where  p  =  density  (mass  per  unit  volume),  w  =  frequency  (radians/sec), 

a  =  side  length,  arid  D  =  the  bending  rigidity  of  the  plate. 

An  important  parameter  in  the  thick  plate  analysis  is  the  shear 

2 

coefficient  a.  In  [40],  Mindlin  suggested  a  value  of  a  to  be  n  / 1 2 . 


In  [16]  and  [37],  a  value  of  5/6  was  used.  When  obtaining  the  results 
of  Table  1,  a  value  of  5/6  was  thus  used  for  u. 
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It  is  seen  in  Table  4.1  that  the  present  results  agree  with  the 
30  elasticity  solution  quite  well  for  practically  all  the  modes  shown 
with  thickness-to-length  ratio  (h/a)  up  to  0.2.  The  maximum  discrep¬ 
ancy  is  2.24%.  The  classical  thin  plate  solution  is  also  shown  in 
order  to  indicate  the  effect  of  shear  deformation  and  rotatory  inertia 
at  the  various  h/a  level. 

Table  4.2  shows  the  first  mode  nondimensional  frequencies  for  the 
same  isotropic  simply  supported  square  thick  plate  but  with  more  values 
of  h/a  and  the  results  are  compared  with  those  given  by  Rao, 
Venkataramana ,  and  Raju  [41]. 


2.2.  Isotropic  Clamped  Square  Thick  Plates 


The  boundary  conditions  were  then  changed  from  clamped  to  simply- 
supported  and  the  results  are  given  in  Table  4.3. 

The  results  by  Rao  et  al.  [41]  was  based  on  a  36  d.o.f.  isotropic 
triangular  plate  element  originated  by  Cowper,  Lindberg,  and  Olson  [25] 
and  Bell  [42].  That  element  was  formulated  on  the  basis  of  Reissner's 
thick  plate  theory.  The  36  d.o.f. 's  include  the  original  18  for 
bending  (see  Figure  4.1)  and  the  additional  six  at  each  of  the  three 
corner  nodes: 


9ip^/ a n 


3  .lip  /3f?,  3\p  /3n 
n  n  n 


where  4*  and  4^  are  defined  as  the  total  slope  due  to  bending  and  shear 
deformation  in  the  fz  and  nz  planes,  respectively.  It  is  noted  that 


Table  4.2  First  mode  nondimensional  frequencies  of  an  isotropic 
simply  supported  square  thick  plate. 


Thickness 

8 

16 

32 

Xdiff- 

Rao  et  al. 

Length 

e  lenient 

e  lenient 

e  lenient 

erence* 

1411 

0.150 


4.002 

3.930 

3.733 

3.410 

3.127 

2.798 


4.000 

3.928 

3.731 

3.408 

3.125 

2.797 


4.000 

3.928 

3.731 

3.407 

3.125 

2.797 


4.000 

3.840 

3.525 

3.184 

2.851 

2.541 


the  present  36  d.o.f.  elements  include  the  18  for  bending  and  the 
additional  six  at  each  of  the  three  corners:  (see  Figure  4.1) 


Ws,  W  .  w 

c  a 


,  w  ,  w5 

en 


nn 


It  is  noted  that  for  the  results  given  in  Table  4.2  and  4.3,  the  shear 
coefficient  a  was  assumed  as  1.0  in  [41]  but  it  was  still  assumed  as 
5/6  in  this  study  in  order  to  be  consistent  with  that  used  in  [16], 
[40]  and  Table  1. 

The  solution  of  [41]  was  obtained  based  on  a  4x4  mesh  with  32 
elements.  It  is  seen  in  both  tables  that  the  present  solutions  has 
converged  at  the  16  element  level.  The  percentage  differences  between 
the  present  results  and  those  of  [41]  increase  as  the  thickness-to- 
length  ratio  increases  and  it  reaches  9.14%  at  h/a  =  0.25  for  the 
simply  supported  plate  and  -8.32%  at  h/a  =  0.25  for  the  clamped  plate. 


ndimensional  frequencies  of  an  isotropic 
ick  plate . 


16  32  Zdiff-  Rao  ec  al. 

element  element  erence  [41] 


14.312 

13.589 

11.817 

9.589 

7.814 

6.246 


13.380 

12.740 

11.156 

9.125 

7.502 

6.039 


13.308 

12.673 

11.104 

9.088 

7.476 

6.022 


13.305 

12.712 

11.258 

9.526 

7.901 

6.523 


*  percentage  difference  as  compared  to  the  results  by  Rao  et  al. 


as  compared  to  the  reaulta  by  the  general ired 
deformation  theory  |44). 


order  plate  theory,  Yang  et  al.  [44]  using  a  generalised  Reissner- 
Mindlin  plate  theory,  and  the  classical  thin  plate  theory.  The  maximum 
discrepancy  between  the  present  results  and  those  using  Reissner-Mindlin 
Theory  is  1.6%. 

The  percentage  discrepancy  for  those  obtained  using  the  thin  plate 
theory  are  given  to  indicate  the  extent  of  the  effect  of  shear 
deformation  and  rotatory  inertia  at  the  various  h/a  ratios. 


2.4.  Three  Ply  (0/90/0)  Orthotropic  Simply-Supported  Square  Plates 
In  Table  4.6,  a  comparison  of  the  lowest  natural  frequency  values 
(«//  E^/ph  for  a  three-ply  (0/90/0)  simply  supported  laminated  square 
plate  using  material  I  in  Table  4.4  is  given.  The  results  for  various 
thickness-to-length  (h/a)  and  beta  values  (e  =  ratio  of  middle-to-outer 
ply  thickness  =  h2  /  hi  =  h2  /  h3;  hi,  h2,  and  h3  are  respectively 
the  thicknesses  of  the  three  plies)  are  compared  with  those  obtained 
by  Bhimaraddi  and  Stevens  [43]  using  a  higher  order  theory,  by  Yang, 
Norris,  and  Stavsky  [44],  and  classical  thin  plate  theory.  It  is  seen 
that  the  present  results  are  in  consistent  agreement  with  those  by 
[44]  based  on  a  generalized  Reissner-Mindlin  shear  deformation  theory 
with  all  the  h/a  and  6  values  considered.  However,  it  is  seen  that  the 
percentage  difference  between  the  present  results  (or  the  solution  of 
[44])  and  those  from  higher  order  theory  [43]  increases  as  h/a  and  B 
increase.  This  discrepancy  reaches  an  extreme  at  h/a  =  0.5  and  e  =  20. 


2.5.  Five  Ply  (0/90/0/90/0)  Orthotropic  Clamped  Square  Plates 

The  natural  frequencies  for  the  first  five  modes  for  a  five  ply 
(0/90/0/90/0)  orthotropic  clamped  square  plate  were  obtained  using 
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three  different  meshes  of  the  present  36  d.o.f.  plate  elements  for  two 
different  thickness/length  ratio  (h/a  =  0.01  and  0.1,  respectively). 

The  results  are  given  in  Table  4.7  along  with  those  by  Craig  and  Dawe 
[45]  for  comparison.  The  nondimensional  frequency  u  is  defined  as 

2  1/2 

a  -  [/-]  (4.11) 

n  yn 

where  is  the  inplane  stiffness  for  the  laminate.  The  present 

analysis  was  based  on  the  same  values  of  the  shear  coefficients 

(oj  =  0.87323  and  a ^  =  0.59139)  as  used  in  [45].  The  same  material 

properties  were  also  used:  E  /E  =  30,  6  /E  =  0.6,  G  /E  =  0.5,  and 

x  y  xy  y  xz  y 

Poisson's  ratio  v  =  0.25.  Each  0  degree  ply  is  two-thirds  the 
xy 

thickness  of  the  90  degree  ply.  In  [45],  a  finite  strip  formulation 
was  developed  and  a  Raleigh  Ritz  method  was  used  to  validate  the  finite 
strip  formulation.  The  results  given  in  [45]  was  based  on  the  Raleigh- 
Ritz  method  using  5  term  series  and  three  finite  strips  using  2  term 
series.  The  results  using  classical  plate  theory  were  also  given  in 
Table  4.7  to  validate  the  accuracy  of  the  present  solution  in  the 
special  case  of  thin  plates. 

The  percentage  differences  between  the  present  solution  and  those 
by  Raleigh  Ritz  method  and  finite  strip  formulation  are  with  2.14%  in 
all  cases. 

2 .6.  Five  Ply  (0/90/0/90/0)  Orthotropic  Simply  Supported  Square  Plates 
The  boundary  conditions  were  then  changed  from  clamped  to  simply- 
supported  and  the  results  are  given  in  Table  4.8  along  with  Craig  and 
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Dawe  [14]  for  comparison.  The  results  given  in  [45]  were  based  on  the 


Raleigh-Ritz  method  using  5  term  series  and  three  finite  strips  using 


2  term  series.  The  results  using  classical  plate  theory  were  also 


given  in  Table  4.8  to  validate  the  accuracy  of  the  present  solution  in 


the  special  case  of  thin  plates. 


The  percentage  differences  between  the  present  solution  and  those 


by  Raleigh  Ritz  method  and  finite  strip  formulation  are  with  2.30%  in 


all  cases. 


2.7.  Anisotropic  (8  =  30  degrees)  Clamped  Square  Plate 


Finally,  an  example  of  anisotropic  (with  fiber  angle  e  =  30 


degrees)  clamped  square  plate  was  analyzed.  The  present  analysis  was 


based  on  the  same  values  of  the  shear  coefficients  (a^  =  =  5/6)  as 


used  in  [45].  The  same  material  properties  were  also  used:  E  /E  = 

x  y 


10,  G  /E  =  0.25,  G  /E  =  0.25,  and  Poisson’s  ratio  v  =  0.3.  The 

a  Jr  Jr  a  Z  Jr  Xy 


nondimensional  frequency  a  is  defined  as 


2  12( 1-v  v  ) 

a_  v  xy  yx_ 

h  Ex 


(4.12) 


The  results  for  the  first  four  mode  nondimensional  natural  frequencies 


■i  for  this  plate  were  obtained  using  various  meshes  of  the  present 
elements  for  two  different  thickness-to-length  ratios  (h/a  =  0.01  and 


0.1).  The  present  results  and  those  obtained  in  [45]  (using  3  finite 
strips  and  a  6-term  series)  are  both  shown  in  Table  4.9  for  comparison 


The  results  show  a  discrepancy  of  about  8.6%  for  the  thickness-to- 


length  ratio  (h/a)  of  0.1  between  the  present  and  the  finite  strip 
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results  by  Craig  and  Oawe.  However,  no  Raleigh  kit/  solution  was  given 
in  [4b J  for  comparison. 

Since  all  the  present  analyses  were  carried  out  using  consistent 
mass  matrices,  it  is  of  interest  to  see  the  results  obtained  using 
lumped  mass  matrix  for  the  same  anisotropic  case  as  shown  in  Table  4.10. 
The  results  show  a  discrepancy  of  about  6.9%  for  the  thickness-to- 
length  ratio  (h/a)  of  0.1  between  the  present  element  and  the  finite 
strip  results  by  Craig  and  Dawe  [45]. 

It  is  noted  that  a  relatively  fine  mesh  of  72  elements  were  used 
and  convergence  of  lumped  mass  results  has  not  quite  yet  been  achieved. 
It  is  expected  with  further  refinement,  the  results  between  the 
consistent  and  lumped  mass  formulation  will  agree. 

2.8.  Anisotropic  Square  Plates 

An  anisotropic  square  plate  with  various  fiber  angles  and  with 
simply-supported  as  well  as  clamped  boundary  conditions  were  analyzed. 
All  the  shear  coefficients  (a^,  a^),  the  material  properties,  and  the 
definition  of  nondimensional  frequency  were  assumed  to  be  the  same  as 
the  previous  example  (4.2.7). 

The  results  for  the  first  four  mode  nondimensional  frequencies  ft 
for  these  plates  were  obtained  using  two  meshes  (4x4,  and  5x5)  of  the 
present  element  for  two  thickness-to-length  ratios  (h/a  =  0.01  and  0.1) 
at  seven  fiber  angles  (e  =  0,  15,  30,  45,  60,  75,  90)  for  the  clamped 
boundary  condition  in  Table  4.11  and  for  the  simply-supported  boundary 


Table  4.10  Nondimensional  frequency  0  for  an  anisotropic 
(33  degrees)  clamped  plate  [lumped  mass]. 


h/a 

Mode 

16  el 

32  el 

50  el 

72  el 

Zdiff- 

H 

Number 

e  lenient 

e lement 

e lement 

e lement 

erence* 

m 

19.896 

20.846 

21.087 

21.171 

0.23 

21.22 

0.01 

26.734 

31.029 

32.263 

32.670 

1 .00 

33.00 

36.482 

42.237 

47.308 

48.999 

2.84 

50.43 

mm 

36.531 

45.898 

49.069 

50.253 

1.87 

51.21 

i 

13.809 

14.653 

14.989 

15.263 

-6.88 

14.28 

0.10 

18.251 

21.137 

22.424 

23.148 

-4.08 

22.24 

20.546 

24.252 

26.099 

27.542 

ao 

c~ \ 

1 

26.64 

KB 

21.969 

26.083 

29.014 

.....  .  j 

30.566 

31  .01 
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percentage  difference  between  present  72  elements  solution 
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SECTION  V 


CONCLUDING  REMARKS 


Three  concluding  remarks  are  drawn  here. 

First,  a  12  d.o.f.  symmetrically  laminated  beam  finite  element, 
the  associated  solution  procedure,  and  a  computer  program  have  been 
developed  for  the  stress  and  vibration  analyses  using  microcomputers. 
The  simplicity  and  efficiency  of  this  development  have  been  evaluated 
through  a  series  of  examples.  The  effect  of  shear  deformation  for  a 
homogeneous  and  isotropic  beam,  the  effect  of  shear  deformation  for  an 
orthotropic  (0/90/0)  laminated  beam,  the  effect  of  shear  deformation 
for  an  anisotropic  laminated  beam  have  all  been  verified  throuyh 
comparison  of  results  with  alternative  existing  solutions.  The  distri 
bution  of  transverse  shearing  stress  through  the  thickness  fur  a  yuasi 
isotropic  laminated  beam  has  also  been  compared.  The  natural  frequen¬ 
cies  of  an  isotropic  homogeneous  plate  and  these  uf  an  anisotropic 
Symmetrical  ly  laminated  cantilever  plate  have  Loti  luipm'd. 

The  program  was  written  in  Basic  1  mynaqe  and  utri<.',  '  i> 

in  the  form  of  half -hand  for  the  static  analysis.  fur  'he  t >  ■  e 
vibration  analysis,  the  Jacobi  eigt-npr  .b  1  •  ,  ;l«i  r  «as  .<  .  la 

expedite'  t  he  erm  put  a  t  ion  ,  the-  pi  uyt  mi  m  !••  ,  i  1  •  •  i  ,  i  ,  i  ,  ’  1  <-r 

A  hard  copy  of  the  plots  1  in  he  attain'd  o.iog  •  >  t.iili'h  >  •  •  n 
dumps."  For  static  unjl  /S  is ,  f  r.e  e'  '<•  ,n  ■"  ■  •  d  r-  t  :  r. 

r.umer  l  c  a  1  data  and  i!  plots  »t  1 1  . '  <  • '  , '  :  or  s  ,  t  :  ■  '  ' 


twisting  angles,  shear  force,  bending  moment,  twisting  moment  along 
the  beam,  ply  stresses  through  the  thickness,  and  the  portions  contri¬ 
buted  by  shear  deformation.  For  free  vibration,  the  program  computes 
the  natural  frequencies,  and  plots  the  mode  shapes. 

Second,  an  18  d.o.f.  symmetrically  laminated  plate  finite  element, 
the  associated  numerical  procedures,  and  the  microcomputer  program  are 
developed  for  the  static,  free  vibration,  and  buckling  analysis  using 
microcomputers.  The  efficiency,  simplicity,  and  practicality  of  the 
development  have  been  evaluated  through  a  series  of  examples. 

The  microcomputer  program  has  i ncorporated  some  highly  efficient 
algorithms  used  in  the  static  and  eigenvalue  analyses.  The  program 
also  gives  a  three-dimensional  plot  of  the  static  deflection  and  mode- 
shapes  of  the  plate. 

It  is  anticipated  that  this  development  could  contribute  to  the 
community  of  structural  analysis  and  designers  who  are  interested  in 
composites  and  u  .e  mi (  rut  mg u t ci s . 

■  t  r  d  ,  if'  i  ,  r  1 1  h  1 1/ r  *  *n ■  i  I  i  .  i  • . f .  i .  t  ?  i  1 1  ,  ,  >  i  , i  ■  ,  l,i  a  r  tli ■  f  it r  n,a  - 

'  '  '  i  ’  ‘  I  •  '  .  ■’  1  V.  a  "  l.o.l.  •  '  i  :  '  a  ii  ,  1  t  .lit  t  <1 

,  '  '  '  I  :  '  ,  i  :  r  ■  !  !  '  .  •  .  •  '  '  '  <  I  •  ■  r  ,  i  '  ’  ■  \  i '  .  \  ■  .  i  i  .  i  I  I  •  1 1 


to  validate  the  accuracy  and  efficiency  of  the  present  element. 
Furthermore,  the  comparisons  served  as  an  indicator  of  the  range  of 
applicability  and  limitations  of  the  present  element. 

In  addition  to  being  used  for  the  purpose  of  application  to  solve 
practical  anisotropic  plate  problem  with  moderate  thickness,  the 
present  formulation  concept  and  procedure  may  be  used  as  a  model  for 
transforming  other  displacement  type  of  isotropic  plate  and  shell 
finite  elements  to  anisotropic  elements  with  the  inclusion  of  the 
effects  of  shear  deformation  and  rotatory  inertia. 
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